Welcome to Bulgaria

Throughout these exercises, you will be working with two different types of satellite imagery for the Kazanlak Valley in central Bulgaria. The images are provided by the JICA and GeoEye Foundation respectively:

• Aster contains the elevation values for the area.
• IKONOS image is a 4-band satellite image covering ca 150 sq km of the Kazanlak Valley with red, green, blue and infrared band, captured in 2001. As it is a fairly large orthophoto, I divided it into two halves with each half at 2 GB and placed them in ScienceData. You can download the full-resolution images manually from public www.sciencedata.dk folder, or directly with with file.download() using these direct links for West and East respectively. In exercise 6 you will play with a slightly reduced East image, which is provided as part of Github data folder.

You will learn to work with these rasters and plot vector data over them.

Task 1: Access raster data values

Raster data can be very big depending on the extent and resolution (grid size). In order to deal with this the raster() and brick() functions are designed to only read in the actual raster values as needed. To show that this is true, you can use the inMemory() function on an object and it will return FALSE if the values are not in memory. If you use the head() function, the raster package will read in only the values needed, not the full set of values. The raster values will be read in by default if you perform spatial analysis operations that require it or you can read in the values from a raster manually with the function getValues().

# Library
___(raster)

# Read in the elevation layer
elevation <- ___("data/Aster.tif")

# Check if the data is in memory
___(___)

# Use head() to peak at the first few records
___(___)

# Use getValues() to read the values into a vector
vals <- ___(elevation)

# Use hist() to create a histogram of the values
___(vals)

rows        2484 
columns     2592 
bands       1 
lower left origin.x        313183.7 
lower left origin.y        4676606 
res.x       30 
res.y       30 
ysign       -1 
oblique.x   0 
oblique.y   0 
driver      GTiff 
projection  +proj=utm +zone=35 +datum=WGS84 +units=m +no_defs 
file        ../data/Aster.tif 
apparent band summary:
  GDType hasNoDataValue NoDataValue blockSize1 blockSize2
1  Int16          FALSE           0          1       2592
apparent band statistics:
   Bmin Bmax     Bmean      Bsd
1 -9999 2347 -3238.349 5106.829
Metadata:
AREA_OR_POINT=Area 
TIFFTAG_DATETIME=2009:04:10 19:05:46 
TIFFTAG_DOCUMENTNAME=created by Program:silcast_cmm 1.09/Ortho at 
TIFFTAG_IMAGEDESCRIPTION=SILC TIFF 
TIFFTAG_RESOLUTIONUNIT=2 (pixels/inch) 
TIFFTAG_SOFTWARE=IDL 6.3, Research Systems, Inc. 
TIFFTAG_XRESOLUTION=100 
TIFFTAG_YRESOLUTION=100 

[1] FALSE

Task 2: Change values: handle missing or bad data values in rasters

There are many situations where you might need to recode raster values. You may want to change the outlier values to NA for example. In the raster package, reclassification is performed with the reclassify() function.

In the elevation raster you’ve worked with the values are meters above sea level and are supposed to range between 0 and 2500. Anything below 0 should be an NA. In this exercise you will assign any values below 0 to NA.

# Plot the elevation layer to see the legend values
___(elevation)

# Set up the matrix
vals <- ___

# Reclassify 
elevation_reclass <- ___(elevation, ___ = ___)

# Plot again and confirm that the legend range is 0 - 2400
___(elevation_reclass)

Task 3: Crop and mask rasters on the basis of other spatial objects

Mask and crop are two similar operations that allow you to limit your raster to a specific area of interest. With mask() you essentially place your area of interest on top of the raster and any raster cells outside of the boundary are assigned NA values. With crop() you limit the extent of your raster to that of your focus area.

The Aster image covers a large area, but we are primarily interested in the areas surveyed by archaeologists, the first of which is the cluster of burial mound points and second, large survey polygons registered by local archaeologists during field survey. In this exercise you will use the extent of mounds and survey units to crop and mask the elevation raster.

On some OSs, raster package does not support sf objects. This should not happen now, but if you encounter difficulty with raster:vector interactions, it helps to know that you can convert the vector to Spatial object with, for example, as(input, "Spatial").

# YOUR CODE

Reading layer `KAZ_mounds' from data source 
  `C:\Users\adela\Documents\RStudio\1_Teaching\cds-spatial\data\KAZ_mounds.shp' 
  using driver `ESRI Shapefile'
Simple feature collection with 773 features and 5 fields
Geometry type: POINT
Dimension:     XY
Bounding box:  xmin: 352481.3 ymin: 4712325 xmax: 371282.4 ymax: 4730029
Projected CRS: WGS 84 / UTM zone 35N

# Read in the survey object
survey <- ___(___)

# Compute the area of the survey
areas <- ___(survey)

# Filter to survey with areas > 30000
survey_big <- ___(survey, ________ > 30000)

Reading layer `KAZ_units' from data source 
  `C:\Users\adela\Documents\RStudio\1_Teaching\cds-spatial\data\KAZ_units.shp' 
  using driver `ESRI Shapefile'
Simple feature collection with 7708 features and 25 fields (with 1 geometry empty)
Geometry type: MULTIPOLYGON
Dimension:     XY
Bounding box:  xmin: 352181.7 ymin: 4712923 xmax: 370892.1 ymax: 4729960
Projected CRS: WGS 84 / UTM zone 35N

# Plot the elevation raster
plot(________)

# Plot the geometry of survey_big
plot(_________(survey_big))

# Mask the elev layer with survey_big and save as elevation_mask
elevation_mask <- ________(elev, mask = survey_big)

# Plot elevation_mask -- this is a raster!
plot(elevation_mask)

Task 4: Reduce the raster resolution (grid cell size) using aggregate()

Rasters, such as ortophotos, terrain models, or aggregated rasters for large area, often come in resolutions far greater than you need. Up till now you have played with fairly small, processed rasters. Now you are getting a taste of the real thing. To reduce the computational load when running analyses, or when developing the right approach, you should use a reduced resolution raster.

The function to reduce resolution in rasters is aggregate() which, as you might guess, aggregates grid cells into larger grid cells using a user-defined function (for example, mean or max). The function used to aggregate the values is determined by the fun argument (the default being mean) and the amount of aggregation is driven by the fact argument (the default being 2).

# Load the IKONOS raster
east <- ________________ 

# Plot the IKONOS raster
___(east)

# Determine the raster resolution
___(east)

# Determine the number of cells
___(east)

# Aggregate the raster
east_small <- ___(east, fact = ___)

# Plot the new east layer
___(east_small)

# Determine the new raster resolution
___(east_small)

# Determine the number of cells in the new raster
___(east_small)

[1] 10 10

[1] 2374596

[1] 100 100

[1] 23940

Task 5: Extract raster values by location

Beyond simply masking and cropping you may want to know the actual cell values at locations of interest. You might, for example, want to know the elevation at each mound location or perhaps the mean elevation within the large survey units. This is where the extract() function comes in handy.

Usefully, and you’ll see this in a later analysis, you can feed extract() a function that will get applied to extracted cells. For example, you can use extract() to extract raster values by point, line, polygon, or neighborhood and with the fun = mean argument it will return an average cell value by neighborhood.

# Extract the elevation values at the mound locations 
________$elev  <- _______::______(elevation, mounds)

# Look at the mounds and extraction results
____(mounds[_______])
hist(____)

Task 6: Raster math with overlay

You will now use the elevation layer and an “prominence” layer. Prominence measures what percentage of surrounding cells are visible from any given location in the raster. A high percentage value thus indicates a prominent point, while a low percentage indicates a low-lying location with poor inter-visibility. Archaeologists assert that mounds sit in high elevations and have a good field of visibility, so let’s see whether they are mostly right :)

What you will do in this exercise is essentially: identify the most prominent locations among the registered archaeological mounds by finding areas that have both a high percentage of prominence and a high elevation. You will use two rasters and you will define an overlay function f to do the raster math with. Remember that raster algebra uses basic Boolean logic to select desirable values. Note that you can only do raster math on two rasters of the same extent, so make sure you align the elevation and prominence accordingly with the crop() function

# Check in on the elevation and read in prominence layer
elevation
prominence <- ___(___)

# Plot prominence
plot(prominence)

# Function f with 2 arguments and the raster math to select specific elevation range and prominence values
f <- function(rast1, rast2) {
    rast1 > _________ & rast1 ________& rast2 ___________
}

# Do the overlay using the above define f function in the `fun` argument. 
# Align the extent of the two rasters with crop() within the overlay() if different

elevation_prom_overlay <- ___(_____(elevation,prominence), ___, fun = ___)


# Plot the result (low elevation and high prominence areas)
___(elevation_prom_overlay)

Task 7 - HOMEWORK: Inspect the most prominent mounds

In the above exercise, we produced an elevation-prominence overlay. Mounds and other sites that sit in this overlay enjoy a strategic position vis-a-vis the rest of the valley. Which ones are they, however, and what are their real prominence values? It is hard to see at the scale of the valley and it would be good to pull the sites out.

Find the mounds that enjoy the most prominent locations as well as those that feature in the elevation-prominence overlay raster. Produce a list of the ID numbers (TRAP_Code) of all the overlay mounds and 25 most prominent mounds and plot them (expressing their prominence somehow) .

# YOUR CODE

Task 8 - HOMEWORK: Sum raster values across polygons - assess the most prominent settlements

Previously you have extracted values from rasters at point locations. More often, you will need to extract values using more complex vector shapes where you will need to decide how to summarize the raster values that fall within the shape.

In this exercise, you will load the layer of ancient settlements for the Kazanlak Valley and calculate and plot the overlapping elevation_prom_overlay values within the settlement polygons. Then you will find out which of the sites were likely defensive in function and contained at least 0.5 ha of strategic high-visibility area.

# Read in the site polygons
sites <- ______("data/KAZ_nuclei.shp")

# Plot the sites
plot(elevation_prom_overlay);plot(__________)


# Create a raster mask from the elevation_prom_overlay using the sites
sites_mask <- _____________

# Crop the raster mask to the sites
sites_strat <- _____________

# Plot and compare the results
plot(sites_mask)
plot(sites_strat)


# Which ones overlap?
defensive <- raster::extract(_______, ______)


# Produce a list of the sites that occupy at least 6 grid cells of most prominent terrain
sites[_______(sapply(defensive, sum)>6),]

Reading layer `KAZ_nuclei' from data source 
  `C:\Users\adela\Documents\RStudio\1_Teaching\cds-spatial\data\KAZ_nuclei.shp' 
  using driver `ESRI Shapefile'
Simple feature collection with 78 features and 57 fields
Geometry type: POLYGON
Dimension:     XY
Bounding box:  xmin: 353091.1 ymin: 4715450 xmax: 370109.9 ymax: 4731864
Projected CRS: WGS 84 / UTM zone 35N

[1] 3227 3231